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Abstract 


We introduce a new family of multivalue and multistage methods based 
on Hermite—Birkhoff interpolation for solving nonlinear Volterra integro- 
differential equations. The proposed methods that have high order and ex- 
tensive stability region, use the approximated values of the first derivative of 
the solution in the m collocation points and the approximated values of the 
solution as well as its first derivative in the r previous steps. Convergence 
order of the new methods is determined and their linear stability is analyzed. 
Efficiency of the methods is shown by some numerical experiments. 
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1 Introduction 
Consider Volterra integro-differential equations (VIDEs) of the form 


y'(t) = g(t, y(t)) +f K(t, T, y(T))dr, tel: (0, T], y(0) = Yo, (1) 


where g € C(S) and the kernel K € C(Q) and satisfies the Lipschitz condi- 
tion with respect to y, with S = {(t,y): te lI,yeER},Q={(t,7y): 0< 
T<t<T, y € R}. It is well known that, under the suitable conditions, (1) 
possesses a unique solution y(t) € C1[0, 7] (see [3]). 
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VIDEs arise in a variety of applications, for example, in viscoelasticity, 
control theory, epidemiology, and population dynamics [10, 18, 26]. 


On the class of discretization methods for the numerical solution of 
VIDEs, Linz [19] proposed algorithms for the numerical solution of (1), which 
consist of linear multistep methods, of the type commonly used for the numer- 
ical solution of the initial value problem for an ordinary differential equation 


y(t) = f(t,y(t)), — y(to) = yo, (2) 


combined with a class of quadrature formulas. In [27], the construction of 
the quadrature rules generated by the backward differentiation formulas was 
discussed in detail and their linear stability properties were analyzed. Also, 
Brunner [2] introduced Runge-Kutta—Nystrém methods, which are based 
on collocation techniques in certain polynomial spline spaces . One of the 
most popular numerical methods for solving this kind of equations is the 
class of collocation methods. In these methods, after discretization of the 
domain, the approximated solution in every subinterval depends on the fixed 
number m of the collocation points. These methods are of convergence order 
m for any choice of collocation parameters; see [5, 3]. Recently, a family 
of multistep collocation methods for (1) has also been introduced, which 
considers interpolation conditions in the previous r step points [7, 8]. Indeed, 
in [12] superimplicit multistep collocation methods (SIMCMs) are used for 
numerical solving VIDEs. In these methods, approximated values of the 
solution in the r previous steps and its first derivative in the m collocation 
points in the current and next subinterval are used, where the r step SIMCMs 
with m collocation points are of convergence order 2m+r. Recently, a method 
based on general linear methods [6] for the numerical solution of (1) has 
been introduced and studied in [20, 21]. Using more of the derivative of the 
approximated solution has been successfully applied to construct methods 
with higher order and extensive region of stability for numerical solution of 
ODEs [17] and nonlinear VIEs [13], specially for stiff problems. 


The purpose of this paper is to construct higher order methods with 
extensive region of stability for solving (1). To do this, we introduce a new 
class of multistep collocation methods, which approximate the solution in 
each subinterval depending on the values of approximated solution and its 
first derivative in the fixed number r of previous time steps, and also the 
values of the first derivative of the approximated solution in the m collocation 
points. These methods will enable us, in deal with stiff equations, to make a 
sensible choice for the steplength of the algorithm. 


This paper is organized as follows: In Section 2, the construction of mul- 
tistep Hermite collocation methods (MHCMs) is described, and in Section 3, 
the convergence orders of these methods are determined. The linear stabil- 
ity is analyzed in Section 4. The paper is closed in Section 5, by showing 
efficiency of the methods by some numerical examples. 
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2 Construction of the methods 


In this section, we describe the construction of the MHCMs for solving (1). 
Let I, = {tn :0=to < ti <--- <ty =T} bea partition of the interval [0, T] 
with constant stepsize h := tn41—tn, n =0,1,...,N—1. In these methods, 
to compute the approximated solution of (1) in the subinterval [t,,,tn+41], we 
use the approximated values of the solution and approximate values of the 
first derivative of the solution in the r previous steps and its first derivative 
in the m collocation points in the subinterval [t,,,t,41]. Denoting the fixed 
collocation parameters by 0 < cy < +--+ < Gm <1 and the collocation points 
by thy =trt+cjh, 7 =1,...,m, the collocation polynomials restricted to 
subinterval [t,,tn+1] are defined by 


r—1 


ery + sh) = yo val 3)Uaa +A Ung thd > xx(s 8)Yn— kp (3) 


k=0 


where s € [0,1], n =r—1,...,N—land U,,; =ul,(tn,;). By differentiating 
from (3) with respect to s, an approximation for y(t) is in the form 


r—-1 
hul, (tn + sh) = yl We + ADL we) Ung thd xk(stn—e (4) 
k=0 


The functions yx(s), ;(s), and vz(s), & =0,1,...,r—l andj =1,2,...,m 
are polynomials of degree 2r+m-—1. By imposing the interpolation conditions 
at points —k, k = 0,...,r —1, in the polynomial (3) and (4) and at points 


cj, j =1,2,...,m, in the polynomial (4), we obtain 
pi(—k) = diz, Yj(-—k) =0, xi (-k) =0, 
1. (C3) = 0, di(cy) = O1j, xi (Cy) = 0, (5) 
yi, (-i)=0, oj(-i)=0, x, (-4) = dix, 

where i,k = 0,1,...,r —1 and l,j = 1,2,...,m. The construction of these 


polynomials is obtained by Hermite—Birkhoff interpolation [24]. Let us as- 
sume the polynomials yx(s), 7;(s) and y;(s) in the form 


2r+m—-1 st 2r+m—-1 2r+m—1 
k 
=D Oils Tale vy? = 2 ar 
(6) 
Now by setting s = c; in (6) and differentiating from these polynomials 
and setting s =c;,, i=1,2,...,mands=-—k, k=0,1,...,r—1, a linear 


system for coefficient of these polynomials is obtained. The coefficient matrix 
Ae ROrt™x@r+m) is in the form 
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1 0 0 0 ies 0 
qk yi (1)? (= 18 (S)) Pree 
im 2 aI **)  G@rtm—i)! 
: ( rt)! ( rt)? ( P41) (—rgay2rtm=1 
1! 2! 3! hia (2r+m—1)! 
0 1 0) 0) suk 0 
2a 42 _y2rtm—2 
Aaie 4 Sr Gem 
5 rey! Cert? (erqayertm=2 
0 1 1! ao ttt (rtm)! 
0 1 ey ct erm 
1! ! (Q2r+m—2)! 
cl C2 c2rtm—2 
0 1 qr a! 0 tt+) rtm) 


and the right-hand vectors are defined by wll € R", k =1,2,...,7, vil) € 
R™, 7 =1,2,...,mas 


[Ky. _ J 9, t#k, uly. _ J 9, ti, 
(u, rath i=k, with. ={ eer 


Their coefficients are obtained by solving the systems 


Aol] = fa 0, Om)”, k=0,1,...,r—- 1, 
Awl = (0,,0,, v7, 3 =1,2,...,m, (7) 
Axl! = [o,,u),0,)7, #=0,1,...,7—1. 


Now, we discuss on the uniqueness of the solution of (7). For this purpose, we 
briefly review some definitions and known theorems about Hermite—Birkhoff 
interpolation (see [15]) 

Let & and n be natural numbers and let 


B=|esls gelocgk, FHV Avan —1 
be a matrix with & rows and n columns having elements 


Ei =0 or 1, 


which are such that 
S- jj =n. 
aj 


We shall also assume that E has no row entirely composed of zeros. Let also 
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G1 <%o< ++ Xp. 

We shall also assume that E has no row entirely composed of zeros. Let also 
Xy< XQ < +++ <r, 

be increasing reals. We also need the set of ordered pairs 

e={(i,j) | ej = 1. 

The reals x; and the “incidence matrix” F describe the interpolation problem 

fO(@) =y! for (i,j) € EB. (8) 


It is appropriate to refer to (8) as a Hermite—Birkhoff interpolation problem, 
which we shall abbreviate to HB-problem. 


Definition 1. We shall say that the HB-problem (8) is poised, provided 
that if 


P(x) € tm~1,P(a;) =0, for all (i,j) € E, 
then P(a) = 0, in the other words, the matrix E is called poised if the 


associated interpolation problem is uniquely solvable for any set of constants 
(3) 


y; , regardless of the choice of the ordered points 2, 72,..., Up. 
Define 
k 1 
m= > G5, M, => mu, j,4=0,1,...,n-1. 
i=l j=0 


Schoenberg [24] showed that a necessary condition for E, to be poised is that 
M,>1+1, 1=0,1,...,n—-1, 
in which these inequalities are called the Polya condition. 


Definition 2. Let the incidence matrix E have k rows. Let f; be the column 
index of the first one that appears in the row 7. Then F is called a pyramid 
matrix if, for each i, €;; = 1 implies e;;, = 1 for f; < j’ <j, and there is some 
value of 1 <i<k such that f; > fo >---> fi and fi < fia < +++ < fr. 


Then Ferguson [15] declared the following theorem for poising the matrix 
FE with respect to the ordering 7, < 42 <-++: < 2x. 


Theorem 1. If & is a pyramid matrix with k rows, satisfying the Polya 
conditions, then E is poised with respect to the ordering x1 < @2 < +--+ < Zp. 
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Now, we show that our interpolation problems (5) have unique solution. 
For these problems, the interpolation points can be considered by 


—r4t1l<-r4+2<---<-1<0<q <@<:::<cm. 
Hence the matrix F with m-+r rows and 2r +m columns can be defined by 
110---0 
110---0 


BH=/110---0 


010---0 
010---0 
Thus, we have 
mo=T, my =m+r, ms = 0, ...,Mn—-1 = 0,7 


Mo =r, M, =2r+m, Mz = 2r +m, oe, Mn—1 = 2r+m=n. 


It can be easily seen that the the matrix F satisfies in the polya condition. 
Also by considering 


fi =9, fo = 0, visg d= 0, fra = 1, oxies Jot = ly 


one can see that the matrix F is a pyramid matrix. Thus by Theorem 1, we 
conclude that F is poised with respect to the ordering —r+1<—-r+2< 
re SC -1L <0 < | < cg < +++ < Gm and the interpolation problems (5) have 
the unique solution. In the other word, the matrix A is nonsingular. 


The exact MHCM is then obtained by imposing the collocation conditions 
for equation (1), that is, the collocation polynomials (3) exactly satisfy (1) 
at the collocation points t,,,;, which leads to the system of m equations in 
the unknowns U,,; in the form 


On = Pai + Oni, (9) 


where 


n—1 1 
Fu = g(tnistin tna) +B Ye | ee ee ee eee 
lo (10) 


o, 4, = nf K(tnji,tn + sh, Un(ty + sh))ds, i=1,...,m. 
0 


Then yn41 = Un(tn41) and hy), = hu), (tr41) are computed by 
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r-1 
Yn+1 = Yvult) )Yn— thw Ung +h Y jrea(d) VO 
"eed (11) 
esr Un— thd vil ng thd xed Uhr 
= j=l k=0 


Also, the discretized MHCM is obtained by using suitable quadrature for- 
mulas for approximating F,,,; and ®,;. The discretized multistep Hermite 
collocation polynomials for approximating y(t, + sh) and y(t, + sh) take 
the forms 


r—1 


Ps (tn + sh) = youl 8)Yn— HALO) Yog +h > xe(o) 8)Yn— ko (12) 


r—1 
hP.( tn + sh) = Yvtte 8)Yn— +a wOo Yj +h dD xe(s) 8)Y,— ko (13) 


where n = r— , N—1 and the unknowns Y,,,; = P/(tn,;) are determined 
by solving the eee system 


By using quadrature formulas with the weights b; and nodes €;, | = 1,..., (41, 
for integrating on [0, 1], and the weights w;, and nodes d;1, 1 = 1,..., uo, for 
integrating on [0,c;], with positive integers io and p,, one can write 


nm-1 ft 
Fag =9(tni, Pr(tn,i)) =f h S- S- bi (tra, ty + ih, Py (ty oh éh)), 
v=0 l=1 
_ Ho 
Gn => wisK (tas tn + dish, Pa(tn + hth). (15) 
l=1 


Substituting from (12) in (15) yields 


m 
Pays =g(tn,i, oa OkYn—k + hy? BigYn gj +h > PikY' n—k) 
k=0 j= k=0 


r-2 M1 


+h x ye wb K (tri, to + Eh, yy (tv + éih)) 


v=0 l=1 


+h S 5 biK Utnisty + ih, Yh thdomy v3 +a Ga k)s 


v=r—1l=1 


HO Al, r-1 


Oni =h wir K (trartn + dish, So FikYn— th Fee ai: k)s 
1=1 k=0 j= b=0 
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where Y,, (t, +sh), v =0,1,...,r—2 are the starting approximated solutions, 
which are obtained by classical collocation method such that 


Ck = pr(c), Bag = Wy), pik = Kel), 
dik = pr(&i), my = 05 (&), Cuk = Xn(E1)) 
Vik = Pr(dit), Bag =Vi(dit),  Pite = Xe(di,1)- 


3 Convergence analysis 


In this section, we study the zero-stability of the method. When the method 
is applied on the equation y’ = 0, the equations in (11) reduce to 


r-1 
Ynt1 = »» Pr(1)Yn—K, 
:=0 


r-1 


hynga = do Pe()yn—e- 
Therefore, in analogy to [8, 22], zero-stability is defined as follows. 


Definition 3. The MHCMs is said to be zero-stable if and only if all roots 
of the polynomials 


p(A) =A" — S- gp(1)Av-*F-3,,  g(A) = AT — S- yi, (LAPP, 
8 k=0 


have module less than or equal to unity and those of modules unity are simple. 


Example 1. In the case, r = 3, m = 2, and c = [ci,1], the method is 
zero-stable for all cy € (0.345, 1]. 


In what follows, we determine the convergence order of the exact and 


discretized MHCMs. We start by deriving local error estimate for the exact 
MHCMs solution P(t) € S$P,,,, (In) for the linear VIDE 


y (t) = a(t)y(t) + g(t) +f K(t,r)y(r)dr, t € [0,7]. (16) 


Using the Peanos theorem [3, 15] for y and differentiating from it, the repre- 
sentation of the local error is obtained as follows. 


Lemma 1. Suppose that p = 2r+m-—1 and that the given functions in (16) 
satisfy a € CP(I), g € C?(I), and K € C?(D) with D := {(t,7): 0<T< 
t < T}. Then for any choice of collocation parameters 0 < cy <+++<¢m <1, 
the exact solution y(t) satisfies 
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vltn + 5h) -¥ eal (tab) +R 4 (8)u! (tg) +2. xe 8)0 tn 
j=l k=0 


+ PRs he: (17) 
and 
hy’ (tn + sh) = y gi.(s)y(tn-e) th wih(s)y’(t 
k=0 j=1 
+Ay xels)u'tt iGice) Tone Ree 8) (18) 
k=0 
with 
Rinrn(s) =f Km,r(s, yt PGs +vh)dy, 
—r+l1 
Kal) =o (3 — 4 =F als)(-k- ye 
k=0 
— (p= 1) 9) vy(s)(-e — 
—h(p—1) y x (8)(—k — ag 
k=0 
and 


_ > 
(2-4 ={",° ts <a t 
) 


Theorem 2. Let e(t) = y(t) — u(t) be the error of the exact MHCM and let 
= 2r+m-—1. Suppose that 


(i) K €C?(DxR), g € C*(d), 
(ii) the starting error is ||é||.0,[0,z,] = O(h”), 
(iii) o(H) < 1, where p denotes the spectral radius and 


w(t ” 

with 

a= Pohl) A> Sth ht, 
A= Shao eo = Lethe eo 
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Then 
llelloo = OCR?" "—*). 


Proof. We point out to the important parts of the proof for MHCMs in the 
case of linear VIDE. The proof can be straightforwardly extended to the case 
of a nonlinear VIDE (1) by using the mean value theorem [3]. Suppose that 
Zn,j = Y (tn,j). By subtracting (3) from (17), representation of the local error 
may be written as 


m r-1 
e(t, + sh) =F) En— rth S-vj(s)E, 7 +h >— xe(SE Z 
j=l k=0 
ethan ia (21) 


with €,_~ = y(tn—k) —Yn—p and Eni = Zp; —U,_;. Differentiating from (21) 
yields 


m T=1 
e' (tn + sh) -¥ ee Ci rth wi(se,,+h>- vi(s) = 
j=l k=0 
+ WR (8). (22) 
Replacing n by J — 1 in (21) and (22), and s = 1, lead to 
EY) = ae + nde, + ne"? + PP Pmrl—-1 (23) 


and 
ne’ — Ae 4 nRev, + nSe'P), + rg eta: (24) 


where A, A, A, and A are given in (20) and 


a 0, 1,m 0,1 1 
S= , , 
(a ). Prana Cain) 
= 0,—1, 0,141 
. ~ ( ( yn) Pm.r -( ) 
w' (1) aa Rag E 
€, = [err41,---,€1]7 € R’, 


1 2 
eV Sane VER. 2 Sanne eR 


(1) = [Yi (1),---,Ym(1)]7, (1) = fi (1), --- din). 


Combining (23) and (24) gives the following matrix equation: 
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EM) ef) (2) 
=H, 7%} + GEM + Qin i_1h”, 25 
(Fe nen, it Qe ai (25) 
where H is given by (19) and 
= hs _ Prrie3 
= (73) , Qn 7 i) , 


The solution of difference equation (25) is (see [16]) 


l-1 
Ep =H HE + SHIGE + HPQmng) (26) 


j=r-1 


On the other hand, setting s = c;, i = 1,2,...,m in (22), by using (4), 
leads to 


Chea =e gO Ben (ea) (27) 


m,rn 


On the other hand, equation (9) for the linear form of VIDE is equivalent to 
n-1 a 
Uni =O(tn,i)Un(tn,i) + (tna) +h > i, K (tn,i, ti + sh)uj(t; + sh)ds (28) 
art 
— nf K(tnistn + sh)un(tr + sh)ds 
0 
and (16) at points t,,; can be written as 
n-1 1 
y'(tn,i) =a(tna)y(tni) + (tna) +h >~ | K (tn i,ti + sh)y(ti + sh)ds (29) 
i=0 "9 


+ nf K (tn istn + sh)y(tn + sh)ds. 
0 


Now by subtracting (28) from (29), we get 


n—-1 


1 
e' (tnt) =a(tn,a)e(tna) +h | Rica eaieae G0) 
eur? 
+ nf K(tnjistn + sh)e(ty + sh)ds. 
0 


Also, by substituting (30) in the (27), we obtain 


n—1 


1 
El, =a(tna)eltn,s) +h Y~ | ee ee eee nee 
1=0 “9 
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+ nf K(tni,tn + sh)e(tn + sh)ds 
0 
_ p2rtm—1 pi (ci). (31) 


m,rn 


By the hypothesis on the starting error and by substituting (21) in (31), we 
have 


r-1 m 
Ee =tGaa).| > nine ta Be), 
k=0 j=l 


r—-1 
+h se Pree =e — inl 
k=0 
n—-1 1 r—1 
+h S- | K (tn i,t + sh) (x pr(s)En—k 
i=0 79 k=0 
m r-1 
+hS° vi(S)Eng thd” xe(S)ER_~ +71" Rm,rn(s) | ds 
j=l k=0 


Ci r—1 
+h | K (tn, tn + 8h) (x yr(8)En—k 
0 k=0 


m r—-1 
+h So vi (s)Eng th Y- xe(S)En_p~ +R” Rmrn(s) | ds 
j=l k=0 


_ pre ae) 
By the hypothesis on the starting error, it follows that 
e(t; + sh) =hPq(s), 1=0,1,...,r—1, (32) 


with ||q||.. < C, independent of h. Hence, we obtain 


n n-1 
(Im — ROMER” =P ST RO +h S71 CED (33) 
1=0 l=F 
thy BOEM + DOE, 
l=r l=r 


where RY eR”, OO, and BY, D® © R™*° are defined as 


7 is K(tnyis ti + sh)aqr(s)ds, 1=0,1,...,r—1, 
(RY), := Se K(tn,i, ti + sh) Rm,r,i(s)ds, l=r,...,2—1, 
Jo K(tn,istn + 8h)Rm,rin(s)ds + a(tn,i)Rm rn (ci), L=n, 
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(B)s. = if Jp, K(tn.as t, + sh)px(s)ds, l=r,...,n—1, 
* K(tnistn + sh)pp(s)ds + a(tni)Px (ci), l= n, 


(co), _ I (tn,i,ti + sh); (s)ds, l=r,...,.n—1, 
: 0 K(tnistn + sh)y;(s)ds + altni)pj (ci), Lan, 


(DW) Vik: aa (tn,isti + sh)xx(s)ds, l=r,...,n—1, 
Jo’ KC tnyirtn + sh)xk(s)ds + a(tni)xn (ci), l=n. 


Substituting (26) in (33), a recurrence formula for E'2) is obtained. Then 
by the same way as described in [9, Theorem 4.2] and by considering the 
starting errors and the fact that p(H) < 1 lead to the estimate 


ED || < Mane, 
and then from (26), a bound for Ex? | is obtained in the form 
EMI] < Mane. 


Using the local error representation and two above inequalities together, com- 
pletes the proof. 


Theorem 3. Suppose that the hypotheses of Theorem 2 hold. If the quadra- 
ture formulas defined in (15) have order 2r + m and 2r +m — 1, respectively, 
then the uniform order of the discretized MHCMs is equal to 2r + m— 1. 

Now, in the table 1, we give a comparison between the methods MCMs 
[8], SIMCMs [12], and the new proposed method (MHCMs) in view of com- 
putational costs with respect to convergence order. 


Table 1: Comparison with multistep and superimplicit multistep collocation 


methods. 
MCMs SIMCMs SIMCMs MHCM 
(type 1) (type 2) 
collocation points m m m m 
dimensional of system m m 2m m 
order m+r—-1 @%wm+r—-1 2m+r—-1 Ar+m-1 


4 Linear stability analysis 


In this section, we analyze the stability properties of the introduced methods 
with respect to the basic test equation (see [5, 3, 4, 23]) 


y'(t) = g(t) + €y(t) + | y(r)dr, t>0, y(0)=%, (34) 
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where €,7 € C. The solution of (34) is stable if Re(A1) < 0 and Re(A2) < 0, 
where Ay,2 = (€ + V/&? + 4n)/2 (see [1]). We observe that, particularly for 
real € and 7, these conditions reduce to € < 0 and 7 < 0. As usual, we look 
for sufficient conditions for the stability of the numerical solution of (34). 


Definition 4. Absolute stability region of the method, FR, is the set of 
all (z := €h,w := nh?) € C x C, such that the numerical solution y, of 
test equation (34) with a constant stepsize h, tends to zero as n + oo. The 
method is Ao-stable if R D R~ x R™ and is A-stable if it is stable for any 
value of (z,w) such that Re(A;) < 0 and Re(A2) < 0. An A-stable method is 
Ao-stable too. 


To state the main results of stability properties of MHCM, let us define 
1 1 1 
OR = | yr(s)ds, Bj = i Vj(s)ds, Yk = 7 Xxk(s)ds, 
0 0 0 


QuK =| pr(s)ds, Vij =| w;(s)ds, Nik =| Vn(s)ds, 
(P(c))ix = pa(ci), (W(e))ig = Yi(ci), (x (€) ix = xa (ci), 


and consider the vectors and matrices 


y(1) =[yo(1),---,r—1(1)]", (1) =[¥1(1),--- Ym (1), 
x(1) =[xo(1),-- + Xr (1), x’(1) =[xo(),---s xr)", 
¢'(1) =[¢0(1),--- ay)’, (1) =[1(1),--- bn)’, 
yl =[Un. +++) Yn—r41 ee ye = 4 ee ere ie 
by Sp tgceenl eal t=ly1..2 4/7 =, 


= [20 J =O]. =O]. eo 


Ey = F g'(1)* | pe = Kea (ees foal (36) 


r—1 0,—-1,1 0;—1,m 0,17 


Theorem 4. Applying the exact MHCM on the test equation (34), leads to 
the following recurrence relation: 


Yn Jin—i ts 
AU, | = R(z,w) | hAU,-1 | +hGy, 
ra v5 


where z := th, w = nh?, 
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-1 
R(z,w) = [Q(z,w)]M(z,w), 


and 
07m 
Gn = O;m ’ 
Gn = Gn—1 
I, 07m 0, 
Q(z, w) = 0, Om I, 
—zp(c) — wT, — z(c) — wP —zx(c) — wA 


FE, F; G, 
M(z,w) = Eo F5 Go 
Mi, Mp M3 


with 
M, =—- zy(c) —-w+ wa, 


Mp =Im — 2(c) — uD + 8, 
M3 = — zx(c) —wA+ wy. 
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Proof. By the notations of this section, we rewrite the first relation in (11) 


in the form 
uns = P(L)TYL? + B(L)PAUn + x(1)P hy’, 
or equivalently 
y= Ey”, + FyhUp-1 + Gihy’),. 
Also, the second relation in (11) can be written in the form 
hy/() = Boy, + FohUy 1 + Gohy', 
Now we apply (10) on the test equation to get 
Uni = Gltni) + €Un(tn,i) + Fr(tni) + On(tn i), 
where it can be written in the matrix form 


Un =Gnté (s(e)u? + (c)hU, + x(c)hy's?) +Fp1+&n, 


with (Gn); = 9(tn,i) and 


n—-1 1 
(Fpa)et= Palisa) = 9h / u(t) + sh) 
i=0 °9 


(37) 


(38) 


(39) 


(40) 
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n—-1 14 fr-l1 m r—1 
= Mh [LS velour n+ YO vil Wis + HT vals) vee | ds 
i=0"%9 \k=0 j=l k=0 


(®,,); = On(tra) = nn ff Un(tn + sh)ds 
0 


r-1 


=nh [ S> yx(s) 5)Yn— Se eee SBT 8) Yh, — k ds. 
0 k=0 j=l j=1 
Thus 
n-1 
Fy-1 = 1h > ay\” + BRU, +-yhy',”, 
1=0 


®,, = nh (ey? 4+ DAU, + Ahy'” ae 
Therefore, we can write 
F,-1 — F,- 2 = = ay”? 1 Se BhUp- 1 + hy’ Lan (42) 


Now for obtaining a recurrence relation, substituting (42) in (41) and multi- 
plying it by h yield 
(In — 20b(c) — wD) AU, + (—29(€) — wy + (2x(e) — wA)hy’,” 
= (Im — 2h(c) — wT + wB) hU,-1 + (—zy(e) — w+ way”, 
t (—zx(c) — wA + why’? + h(Gn — Gn-1): 


(43) 


Hence the relations (38), (39), and (43) complete the proof. 


Let us define 
M1 a Hi H1 
Ge => — dbipe(&), By = >_bvs(E), Fe = d_ bixe(G), 
1=1 I=1 1=1 


Ho Ho 


a= Wil Pr ( dit) te 2 wi; (d 4,1) Bu => WilXk( dit) 


Theorem 5. The discretized MHCM, applied to the test equation (34), leads 
to the following recurrence relation: 
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y ss Yn-1 
hY, | = R(z,w) | hYn-1 | +hGa, 
ye yo 


where z := th, w = nh?, 


R(z,w) = (Q(z, w)|-!M(z, w), 


and 
I, 0; 0,7, 
Q(z, w) = 0,7 _ Om = I, _ ; 
—zp(c) — wOI,, — z(c) — wP —zx(c) — wA 
_ E, Fi Gi 
(ew) = | Ea Fa Gs 
M, My M3 
with 


M, =— zp(c) — w+ wa, 
Mz =In — zap(c) — uD + wB, 
M3 =—- zx(c) — wA + wy. 


Proof. It is similar to the proof of Theorem 4. 


The R(z, w) is called the stability matrix of the method. Now, the method 
is stable if p(R(z,w) < 1 and the stability region of the method is R = 
{(z,w) € Cx C: p(R(z,w) < 1}. Here, the term G,, does not influence 
stability. The stability function of the method with respect to (34) is defined 
as 


p(z,w;A) = det(Alor4m — R(z, w)). (44) 


To investigate the stability properties of the exact MHCM, it is more con- 
venient to work with the polynomial obtained by multiplying the stability 
function (44) by its denominator. The resulting polynomial will be denoted 
by the same symbol p(z, w; A). This polynomial takes the form 


2r+m 
p(z,w;A) = S- pi(z,w)’, (45) 
i=0 
where p;(z,w), i = 0,1,...,27 + m are polynomials of degree less than or 


equal to m. Denoting the roots of the polynomial p(z,w;A) by Ai, 2,-.-, 
A2r+m, the absolute stability region of the method is then defined by 
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Figure 1: Stability region with r = 2, m= 2, c= (s, 1]. 
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Figure 2: Stability region with r = 3, m= 2, c= [52, 1]. 


R={(z,w)€eC xC : |A,(z,w)| <1, +=1,2,...,2r+ m}. 


We did not find Ag-stable methods within this class, but wide stability 
regions exist. For example, in the cases r = 2, m = 2 with collocation 
parameters c = [4 1], and r = 3, m = 2 with collocation parameters c = 
a 1], the regions of stability are unbounded. These regions are given in 


Figures 1 and 2, respectively. 


5 Numerical results 


In this section, to check the numerical performance of the method, we have 
considered a variety of problems. Here, the starting values y;,...,y,—1 are 
obtained by a one step MHCMs of the same order of the present method. 
In practice, we need quadrature rules to obtain the numerical solutions. 
For this propose, we have to apply the rules that preserve order of the main 
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Figure 3: Stability region with r = 3, m=2, c= (2, 1). 


method. A suitable choice is the 3-times Romberg quadrature formula with 
the Simpson rule [25]. In this rule that is of order 8, the nodes and weights 
for integration in the interval [0,h] are 


and 


ih 
G=F 20,1,..458 


1 


>= =e0l 


217 1024 352 1024 436 1024 352 1024 217]. 


In what follows, we describe details of the implemented methods: 


Method 1: MHCM of convergence order 5 with r = 2, m = 2, and 
collocation parameters c, = 2 and cp = 1. 


Method 2: MHCM of convergence order 7 with r = 3, m = 2, and 
collocation parameters cj; = 3 and cg = 1, with bounded stability 
region. (see Figure 3) 


Method 3: MHCM of convergence order 7 with r = 3, m = 2, and 
collocation parameters cy = ae and cg = 1, with unbounded stability 
region. 


Method 4: Multistep collocation method [8] of convergence order 4 
with r = 3, m = 2, and collocation parameters c, = 5 and co = 1. 


Computational experiments are done by applying the Methods 1-4 on the 
following problems: 


I. 


The linear VIDE 


t 
y(t) =1 + 2t — y(t) + i r(1+ Ar)eT My (r)dr, t € (0, 2], 
0 
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Table 2: The results for problem I. 


kK [4 5 6 7 8 9 
Method 1 cd [2.74 420 569 718 868 10.19 
p(h) 4.86 4.93 4.97 4.98 4.99 
Method 3 cd | 4.02 6.06 814 10.23 12.33 14.44 
p(h) 6.77 6.90 6.96 6.98 6.99 
Method 4 cd | 2.94 4.07 5.22 641 760 881 
p(h) 3.71 4.84 3.92 3.96 3.98 


y(0) =1, 
with the exact solution y(t) =e’ . 


II. The nonlinear VIDE 


Phi 1 | 1, 242, f° dr 
y()}=—t (+2)? | mn ) ‘| iLife (0, 4], 
y(0) =1, 


with the exact solution y(t) = TH 


III. The stiff nonlinear VIDE 


j r A, oA, — 
y(t) = 5 ult) F1+ 5t+ ste A} tre¥ dr, +t € [0,4], 
0 


with the exact solution y(t) = t. 


We have implemented the methods with a fixed stepsize h = x, with 
several integer values of k. In the following tables, the maximal end point 
error is written as 10~°¢, where cd is the number of correct significant digits. 
Also, a numerical estimation of the order of convergence of the methods 
is computed by the formula p(h) = log, (GY), where e(h) is the maximal 


absolute end point error. 


The results in Tables 2 and 3 confirm the proved convergence order. In 
Table 4, we show the effect of linear stability properties of the methods in 
solving stiff problems. For Method 2, the interval of absolute stability on real 
line is bounded. When z = Ah or w = Ah? lies out of this interval, increas- 
ing of absolute error is evidently seen while for Method 3, with unbounded 
stability region, better results are obtained. 
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Table 3: The results for problem II. 

k 5 6 7 8 9 

Method 1 cd 6.45 7.85 9.31 10.79 = 12.28 
p(h) 4.67 4.84 492 4.96 

Method 3 cd 747 9.33 11.31 13.36 15.43 
p(h) 6.17 6.57 6.78 6.89 

Method 4 cd | 461 567 6.77 7.92 9.06 
p(h) 347 3.65 3.82 3.91 


Table 4: Comparison of the methods (absolute errors of the methods for problem III 
with \ = 500). 


Method 2 Method 3 

t k=7 k=8 k=7 k=8 
0.25 | 1.44E-20 1.28E-22 | 3.25E-20 1.38E-22 
0.5 5.69E-19 8.35E-22 | 2.14E-19 8.58E-22 
0.75 | 2.20E-17  1.89E-21 | 4.79E-19 1.90E-21 
1 6.21E-16 2.93E-21 | 7.34E-19 2.90E-21 
1.5 4.86E-13 7.34E-21 | 1.82E-18 7.34E-21 
2 4.24E-10 1.71E-20 | 4.28E-18  1.72E-20 
3 3.47E-04  3.09E-20 | 7.72E-18 3.09E-20 
3.5 3.13E-01  3.57E-20 | 9.03E-18  3.61E-20 
4 2.73E+02 4.13E-20 | 1.03E-17 4.13E-20 


6 Conclusion 


The introduced methods for VIDEs, based on Hermite—Birkhoff interpolation, 
use of the approximated values of the solution in the m collocation points 
and the approximated values of the solution as well as its first derivative in 
the r previous steps. Using of this technique caused to get methods of higher 
orders and also with extensive stability regions. As we showed in Table 4, 
the extensive stability region of the method let us to make a sensible choice 
for the steplength of the algorithm for solving stiff equations. 
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